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Dynamical mean field study of the Mott transition in the half-filled Hubbard model 

on a triangular lattice 
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We employ dynamical mean field theory (DMFT) with a Quantum Monte Carlo (QMC) atomic 
solver to investigate the finite temperature Mott transition in the Hubbard model with the near- 
est neighbor hopping on a triangular lattice at half-filling. We estimate the value of the critical 
interaction to be U c = 12.0 ± 0.5 in units of the hopping amplitude t through the evolution of 
the magnetic moment, spectral function, internal energy and specific heat as the interaction U and 
temperature T are varied. This work also presents a comparison between DMFT and finite size 
determinant Quantum Monte Carlo (DQMC) and a discussion of the advantages and limitations of 
both methods. 
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I. MOTIVATION 

Systems with triangular lattice structure have been 
a source of attention mostly due to the frustration ef- 
fects resulting from their non-bipartite structure. As 
a result of the competition between the frustration and 
strong electron correlations, these systems exhibit a wide 
range of exotic phases. Recent studies of the metal- 
insulator transition, superconductivity, and antiferro- 
magnetism in the organic compounds k-(BEDT-TTF)2X 
with X as an anionic, discovery of superconductivity in 
Na x Co02.yH 2 0^ and the recent discovery of the Mott 
transitions in 0.33 monolayers of Sn on Ge(lll) at 30 
K^ as a few examples, underline the importance of these 
systems and their physics. 

Theoretical work has been dedicated to the discovery 
of the magnetically ordered phases in the ground state of 
the Hubbard and t — J models on triangular lattices as a 
function of the on-site electron-electron Coulomb inter- 
action U at different fillings jJi&LiuLiJl Unlike the square 
lattice at half-filling, which is a Mott insulator with an- 
tiferromagnetic order at arbitrarily small values of U/t, 
the ground state of a triangular lattice has a variety of 
magnetically ordered and disordered phases. This is due 
to the lack of perfect nesting in the non- interacting Fermi 
surface of a triangular lattice at half-filling. 

For triangular lattices, according to the Hartree-Fock 
calculations of Krishnamurthy and co-workers p£ the 
Mott transition occurs from a paramagnetic metal to 
a paramagnetic insulator at half-filling for values of 
U larger than the band width W = 9t. A variety 
of the physical properties of the triangular lattices in- 
cluding tendencies towards superconductivity within the 
small to intermediate U regime have been studied at 
finite temperature using correlated electron approaches 
which go beyond mean field theory, such as the fluc- 
tuation exchange approximation (FLEX)ii or one-loop 
renormalization-groupJ^ However, when U/W > 1, a 
more powerful cluster solving technique such as the 
Quantum Monte Carlo (QMC) is required to accurately 
describe the phase transition. Finite size lattice determi- 



nant Quantum Monte Carlo (DQMC) has already been 
employed for triangular and kagome latticesiiii4 Unfor- 
tunately, the finite size lattice QMC method incurs sign 
problems at both low temperatures and away from half- 
filling. In addition, due to the finite size nature of the 
problem, the non-local correlations are always overesti- 
mated and the system consisting of a finite size cluster 
becomes an insulator as soon as the correlation length 
reaches the size of the lattice at low enough T. Thus, 
the gaps or pseudogaps are also overestimated and are 
not guaranteed to persist as one reaches the thermody- 
namic limit. 15 However, finite size effects can be reduced 
by size scaling and in general become less important as 
U increases. 

In this work, we employ the dynamical mean-field the- 
ory (DMFT) 16 ! 17 : 18 to study the Mott transition for the 
Hubbard model on a triangular lattice at half-filling. 
DMFT describes systems in the thermodynamic limit 
and at the same time suppresses the physics due to the 
non-local correlations through coarse-graining over all 
momentum degrees of freedom. Therefore, susceptibil- 
ities can diverge, as in the thermodynamic limit but the 
Mott gap in the strong-coupling regime is always smaller 
than its counterpart given by a finite size lattice ap- 
proach. The combination of the DMFT and QMC also 
does not incur sign problems at either low temperatures 
or any fillings which makes the low temperature physics 
more accessible and reliable. 

DMFT is a more suitable approximation for the tri- 
angular lattice systems as opposed to those with square 
lattice structure. For instance, DMFT is incapable of de- 
scribing the metal-insulator transition on the square lat- 
tice in the weak coupling regime. This is a consequence of 
the suppression of the short-range antiferromagnetic cor- 
relations responsible for this transition at small U by the 
single site nature of DMFT. I 5 . However, in triangular lat- 
tices, due to the frustration of the lattice structure, even 
in the presence of non-local correlations, the short-range 
antiferromagnetic correlations are greatly suppressed. As 
a result, heavy quasiparticles form at the Fermi energy 
reminiscent of the DMFT physics^ 



DMFT has already been employed for studying trian- 
gular lattices using the exact diagonalization and Lanc- 
zos techniques &L In this paper, we present for the first 
time, the combination of DMFT and QMC for solving 
the Hubbard model at half-filling on a triangular lattice. 
By computing the density of states at different values of 
interaction U and temperature T, we estimate the value 
of the critical interaction for the Mott transition to be 
U c /t = 12.0 ± 0.5, less than the predicted value of 15 
in Ref. l2fl but, as expected, larger than the mean field 
result. In addition, the Mott transition features are also 
manifest in the evolution of the magnetic moment, to- 
tal internal energy and the specific heat as functions of 
on-site interaction U and temperature T. 



II. MODEL AND FORMALISM 

We consider the repulsive Hubbard model on a trian- 
gular lattice, 



with t the hopping amplitude, fj, the chemical potential 
and U the on-site Coulomb interaction between the elec- 
trons of opposite spins residing on lattice site i. Note that 
unlike the square lattice, on a triangular lattice, due to 
the lack of particle-hole symmetry, \x = U/2 does not 
yield half-filling. Also, this model exhibits quite differ- 
ent physics under the change of the sign of t, especially 
away from half-filling. The dispersion in the U = limit 
including the nearest neighbor hopping only is, 

k \/3 
e k = -2t{cos(k x ) + cos(-^ + — k y ) + 

cos(^ - ^k y) ] . (2) 

The resulting noninteracting density of states is markedly 
asymmetric about E = 0, extending from E = — 6t to 
E = +3t, with a van Hove singularity at E = t. 

The DMFT approach to tight binding Hamiltonians 
starts by coarse-graining the Green's function, 
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(3) 



where u) n = (2n + 1)ttT is the Matsubara fermionic fre- 
quency and S cr (iw n ) the local self-energy for spin a. Note 
that the self-energy T, a (iuj n ) in DMFT is local (momen- 
tum independent). In the first DMFT cycle iteration, we 



chose T, a (iLu r , 
computed, 



0. The bath Green's function is then 



Q a {Wn) = (G^fan) + Y*cr{iuj n )y 



(4) 



Qu{iijJ n ) is used as an input to the QMC solver. 
Our DMFT-QMC solver is based on the Hirsch-Fye 
algorithm2ii22 in which the interacting two-body term 
in the Hubbard model is decomposed into two one-body 
density terms coupled to a Hubbard-Stratonovich field. 
Through the stochastic averaging over the most probable 
configurations of these Ising like fields, one calculates the 
fully interacting Green's function on a single site of a lat- 
tice G a {iuj n ). The local self-energy E CT (iw„) is then cal- 
culated from the full and bath Green's functions G a (iujn) 
and Ga{iuj n ) respectively, 



S CT (iu> n ) = Q„ 1 (iw„) - G a 1 (wi>„) . 



(5) 



This self-energy is again inserted in Eq|3 and the self- 
consistent process continues until convergence. 

The final product of the DMFT-QMC cycle is the full 
Green's function in imaginary time G(t) which is used 
to compute the density of states (DOS) N(u>), 
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N(uj) must be calculated by inverting Eq|Sl which is an 
ill-posed problem due to the statistical errors in G(t) 
from the QMC simulations. This task in this paper 
is done by using the maximum entropy (MaxEnt) an- 
alytical continuation technique developed by Jarrell and 
Gubernatis. 2 - This technique is based on Bayesian in- 
ference in which an entropy function with respect to an 
original default model is defined as a function of N(u>). 
The best N(ui) is the one that maximizes this entropy 
function for a given G(t). 

We also present results for the double occupation, 



< D >=< n^nii > , 



(7) 



which is suppressed by the on-site repulsion. The double 
occupation is related to the magnetic moment, 



< n iT > + < mi > -2 < D > 



(8) 



The total internal energy can also be computed in this 
formalism in terms of the full Green's function and double 
occupation: 

E mt = ^ 5> fe G CT (£;, zc J „)]e^" 0+ + UD , (9) 

nha 

with D defined as in Eq[7| Eq|5] is equivalent to 
the Galitskii-Migdal expression for the total internal 
energy^ The specific heat at constant volume is 
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,dE, 



dT ' v 



(10) 



As seen later in this article, the variation of the total 
internal energy and specific heat as a function of U and 
T will exhibit features reflecting the Mott metal-insulator 
phase transition. 



III. RESULTS 

Before presenting our data for the density of states, we 
comment on some of the details of the maximum entropy 
procedure, since careful characterization of the data is 
essential to obtain reliable results. In this work, the 
DMFT-QMC cycle as described in section [H] converges 
at a desired tolerance which is about few percent for the 
time Green's function G(t) . The chemical potential is 
also adjusted simultaneously to keep the total occupa- 
tion at half-filling. For the QMC runs, we adopt the 
criterion U(At) 2 < 0.05 with At = 0/L for finding L, 
the number of the time slices. L can be small at high 
temperatures, however, even for temperatures as high as 
the bandwidth, we choose the lower limit of L = 40. The 
number of the lattice sweeps for the QMC runs varies be- 
tween N = 2 x 10 5 to 1 x 10 6 depending on the values of 
U and T and the reliability of the MaxEnt results. The 
outputs of the cycle are the fully interacting time Green's 
function, the magnetic moment on a single site and total 
internal energy for the temperature and interaction at 
which the QMC cycle has been running. The full G(t) is 
employed to calculate the DOS through EqEl using the 
MaxEnt technique^ 

Care must be taken when using MaxEnt as the results 
can strongly depend on the quality of the QMC data. 
MaxEnt, as described in Ref.l23l requires statistically in- 
dependent and Gaussianly distributed measurements of 
G(t) for every time slice r from QMC. Both these two 
properties are often lacking in the G(t) data initially ob- 
tained from the QMC output because of the presence of 
correlations between the measurements. The solution is a 
careful rebinning of the G(t) data which eliminates these 
correlations and yields Gaussianly distributed data. Cor- 
relations between different time slices are also eliminated 
in the MaxEnt code by diagonalizing a covariance matrix 
coupling different time slices together £& 

The other concern is a proper choice for the MaxEnt 
default model. In general, for good QMC data, the choice 
of the default model should not change the qualitative 
features of the spectral function significantly. Among 
the available default models, the Gaussian and second 
order perturbation models have been most commonly 
utilized. In our work, we adopt the latter and perform 
the so called annealing method in which we start with 
the second-order perturbation theory default model at 
a high temperature close to the band width. We then 
use the MaxEnt output for that temperature as the de- 
fault model for a slightly lower temperature and continue 
until we arrive at the lowest desired T. Finally, as intro- 
duced in Ref. |23, there are different methods of doing 
the MaxEnt including the historic, classic and Bryan. 




FIG. 1: (color on-line) The evolution of the density of states 
at half-filling for the triangular lattice as a function of U . For 
clarity, the base lines have been shifted in steps of 0.1 for 
different values of U. The ui = line determines the location 
of the Fermi energy. Three behaviors are visible: At weak 
coupling the density of states has a finite value at uj = 0. 
At large coupling the density of states vanishes, signaling the 
opening of the Mott insulating gap. At intermediate coupling 
the density of states exhibits a Kondo resonance at u) = 0. 



We utilized both classic and Bryan techniques in our cal- 
culations and found that both these two methods give 
qualitatively similar results. 

In FigQ] we present the evolution of the DOS as a func- 
tion of U at low temperatures. At U = 0, the DOS shows 
a Van Hove singularity near u) = t. As the interaction 
is gradually turned on, there is an overall broadening of 
the density of states, including a smearing out of the Van 
Hove singularity. The DOS at U = 8£ and T = 0.2t ex- 
hibits the onset of the Abrikosov-Suhl (Kondo) resonance 
in which there appears a sharp quasi-particle peak near 
the Fermi energy with two shoulders, the upper and lower 
Hubbard bands, around it. For U = lOt and T — 0.2t, 
the peak lines up with the Fermi energy and shoulders are 
approximately located at ±5£ « ±Z7/2. At U = 12t, and 
T = 0.25£, we clearly observe a different trend in which 
the Kondo resonance is absent and the DOS has a pro- 
nounced minimum at the Fermi surface, indicating that 
the system has been driven to the metal-insulator transi- 
tion regime. For U = 14£, and T = 0.5t, the (Mott) gap 
is well established and the system becomes an insulator. 

To better address the role played by temperature in the 
metal-insulator transition, in Fig|2]the DOS for three val- 
ues of U = 10£, 12t and lit and different temperatures 
have been plotted. At U — lOt, as T decreases, the 
DOS evolves towards a sharp Abrikosov-Suhl resonant 
quasi-particle peak. However, U — Y2t exhibits a com- 
pletely different trend. As the temperature drops, the 
DOS evolves towards the opening of a gap. At T = Q.25t, 
we noticed that the MaxEnt results are of less reliabil- 
ity due to the correlations in the QMC data which could 
not be totally removed up to 10 6 QMC lattice sweeps by 
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FIG. 2: (color on-line) Panel (a): the DOS for U = lOt at 
half-filling and different temperatures. As the temperature is 
lowered a sharp Kondo resonance develops at u> = 0. The 
system is metallic. Panel (b): the same results as for panel 
(a) for U = 12t. The system is on the verge of insulating 
behavior. Panel (c): U = 14t. A Mott gap has opened at 
T = 0.5i. 




1 

0.8 




" 


0.6 
0.4 


I : 


0.2 



- 


2 4 6 8 10 12 14 


U 




•-•U=4t 




G-O U=4t (DQMC) 




»-«U=6t 




♦— »U=8t 




O-O U=8t (DQMC) 




*— »U=10t 




<-^ U=12t 




U=14t 





10 



FIG. 3: (color on-line) Variation of the magnetic moment 
with temperature for different values of interaction U at half- 
filling. < m'i > drops down close to its minimum 0.50 at large 
temperatures (cut off in this figure). The comparison between 
the DQMC and DMFT at U = At and 8t is also presented. 
The inset shows how T ma x , the temperature for the maximum 



< m z >, drops as U increases. 



rcbinning the data. (It is, of course, typical of Monte 
Carlo simulations that statistical fluctuations are largest 
near phase transitions) This is additional evidence that 
U = 12£ is close to the critical value for the opening of 
a gap but we were not able to determine conclusively 
whether the gap exists at U = 12£ if one can reach low 
enough temperatures. On the other hand, by the time 
U = 14i (at T = 0.5£), the MaxEnt results show a pro- 
nounced Mott gap in the DOS. Therefore, it is clear that 
the critical value of U c must lie very close to U = 12t, 
leading to our estimate U — 12t ± 0.5i. This value is 
smaller than the value U c « 15t obtained in Ref. |2(J us- 
ing the DMFT with exact diagonalization and Lanczos 
techniques as the impurity solvers, and comparable to 
U c « 12t in Ref. 7 with exact diagonalization on 12 site 
lattices. At both values of U — 12t and 14£, we also ob- 
serve the appearance of a feature on the right hand side 
of u) = 0. This feature is likely a remnant of the peak in 
the non-interacting density of states- it disappears as U 
is increased to larger values and is completely absent for 
U = 16£. 

Another signature of the metal-insulator (Mott) tran- 
sition is manifest in the behavior of the magnetic moment 
defined in Eq|S]as a function of temperature T and inter- 
action U. As shown in Fig|3 for small to intermediate 
values of U, < m 2 z > reaches a maximum at a value of 
Tmax which is characteristic of a Fermi-liquid dominated 
by spin fluctuations at low temperaturesii£*2£i2£ At low 
temperatures with an entropy of jT per particle, a Fermi 
liquid gains free energy upon heating by trying to localize 
the electrons in order to take advantage of a larger spin 
entropyi2Si2i The inset in Fig[3|shows how T max decreases 
with U. However, at the onset of the Mott transition, 
< m 2 z > becomes a monotonically decreasing function of 



T as seen in Fig. [21 for U = 12t and lAt. In a localized 
Mott phase, the spin entropy at T ss is finite and the 
system cannot gain free energy by further localizing the 
electrons upon heating. Comparison has also been made 
between DQMC on a 6 x 6 finite size lattice and DMFT 
in FigOHfor values of U = At and 8i. The agreement be- 
tween the two approaches is very close at U — At down 
to T — 0.125£. However, for U = 8t, at low tempera- 
tures, the DQMC approach exhibits a rather flat region 
for < m 2 z > below T fa It in place of the downward curve 
in DMFT. This flat feature is similar to what is observed 
at larger values of U in DMFT as a precursor to the in- 
sulating phase. This is typical of any finite size lattice 
approach such as DQMC in which at low enough temper- 
atures, once the correlation length exceeds the size of the 
lattice, the system turns insulating. In the DMFT, on the 
other hand, the system is in the thermodynamic limit and 
therefore, the insulating phase features do not appear ex- 
cept at large enough interaction U values. Nevertheless, 
DQMC has the advantage of including the non-local cor- 
relations which are ignored in DMFT. These correlations 
work in favor of the gap formation. Thus, DMFT over- 
estimates the value of U for the formation of the gap in 
the Mott transition. 

The evolution of the total internal energy and specific 
heat defined in Eq 9 and Eq 10 respectively also exhibits 
signatures of a metal-insulator phase transition at large 
interaction U values and low temperatures T. In Fig01 
panel (a), the total internal energy has been plotted as 
a function of temperature T for different values of inter- 
action U. We have employed a polynomial fitting proce- 
dure developed by Duffy and Moreo^ and obtained the 
curves presented in panel (b). By analytically differenti- 
ating the polynomial fits with respect to temperature, we 
arrive at results for the specific heat plotted in panel (c). 
For values of U = At to lOi, we observe a two-peak struc- 
ture in the Cy consistent with the results in Ref. I27l28l 
for the infinite dimensional and 2D hypercubic Hubbard 
models respectively. The peak at low temperatures cor- 
responds to the local spin fluctuations analogous to the 
formation of the Kondo singlets in the Kondo model. It 
is also a signature for the validity of the Fermi liquid pic- 
ture as it indicates the existence of a conducting electron 
bath around the magnetic impurity^ The upper peak 
at high temperatures, however, corresponds to on-site 
charge fluctuations. Separation between these two peaks 
occurs at U = At and develops up to U = 8£ resulting in 
a sharp lower peak. 

As prescribed in Ref. |23, the lower and higher temper- 
ature segments of Ei nt (T) have been fit with two different 
polynomials which meet around T fa It. This prescrip- 
tion explains the rather sharp feature in the Cy curve at 
U = At in Fig|H panel (c), around T fa 1.25t. However, 
the necessity of using two polynomials originates in the 
existence of two separate temperature regions, one dom- 
inated by spin and the other by charge fluctuations. The 
lower peak dwindles at U = lOt as a sign of approaching 
the insulating phased It is also striking that, similar to 



the result of Ref. 1271 within the Fermi liquid range up 
to U = lOt, all Cy curves more or less intersect at the 
same temperature close to T = 2. It. The picture how- 
ever is different at U = lAt when the insulating phase 
sets in. The lower peak completely vanishes and the Cy 
curve goes to zero at a finite value of T. This is consis- 
tent with the behavior of Cy in the presence of a gap A 
going as e~ ' T at low temperatures, and having a broad 
peak at high temperatures resulting from charge fluctu- 
ations in the upper Hubbard bandit! For U = 12£, E int 
results were too noisy to be conclusive and therefore are 
not presented. This was apparently due to U = 12t being 
close to the critical U c which made it difficult to achieve 
convergence in the simulation. 

Lastly, Fig[S] presents a comparison between DMFT 
and DQMC for E lnt {T) at U = At and 8t. Similar to 
the trend observed for < m 2 z > in Fig|3| at high enough 
T, the agreement betweeen the two approaches is clearly 
manifest in the results. For U = At at low T, DQMC 
exhibits an insulating phase feature in the rather flat tail 
of the Ei nt (t) curve near T — 0, typical of a finite size 
approach as discussed earlier. For U = 8t, however, the 
effect of non-local correlations becomes more important 
and at the same time, finite size effects diminish. Hence, 
DQMC features may be more realistic in showing insu- 
lating phase signatures. Nevertheless, U = 8t is still 
below the bandwidth W — 9t and we do not expect to 
have entered the insulating phase as yet. Unfortunately, 
mostly due to the sign problem in DQMC, attempts in 
going to lower temperatures and higher U values have 
not yet been successful. Thus, it is not clear to us which 
approach more plausibly describes the physics for the in- 
termediate value of U = 8t. 

In order to bring the DQMC and DMFT treatments 
into full agreement, it would be necessary to extrapolate 
both to a zero value of the discretization size of the in- 
verse temperature (3 and also to extrapolate DQMC and 
DMFT to infinite spatial cluster size and full momentum 
resolution (the dynamical cluster approximation (DCA)), 
respectively. 



IV. CONCLUSIONS 

In summary, we have performed the first implementa- 
tion of the DMFT using QMC as the solver for the half- 
filled triangular lattice Hubbard model. Our approach 
confirms the previous predictions that for large enough 
U/W, there is a metal-insulator (Mott) phase transition 
from a paramagnetic metal to a paramagnetic insulator 
at finite temperature. We demonstrate this phase tran- 
sition by presenting the behavior of the DOS, magnetic 
moment on a single site, total internal energy and spe- 
cific heat as functions of the interaction U and temper- 
ature T. Our results suggest a critical value of interac- 
tion U c /t = 12.0 ± 0.5. This is consistent with the value 
U c sa 12t from exact diagonalization on 12 site lattices^ 
and U c fa 15£ from DMFT with exact diagonalization 
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teraction U including the error bars. Panel (b): Polynomial 
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U taken by analytical differentiation of polynomials in panel 
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and Lanczos impurity solversiSi. We have also investi- 
gated the intermediate interacting regime which is char- 
acterized by the presence of an Abrikosov-Suhl (Kondo) 
resonance. 

DMFT is a more suitable approximation for studying 
the Mott transition on a triangular lattice compared to a 
square lattice. The reason is that in a triangular lattice, 
due to the frustration of the lattice structure, the short- 
range antiferromagnetic correlations are suppressed and 
do not play a major role in the metal-insulator transi- 
tion at small U*& Hence, neglecting these correlations, as 
done by DMFT, should still give us a qualitatively correct 
physical picture at least in the large U limit. The combi- 
nation of the DMFT and QMC also does not encounter 
the difficulties related to the sign problem at low tem- 
peratures and away from half-filling and also finite size 
effects associated with DQMC at low T. Thus, the low 
temperature physics is better described by the DMFT 
due to being in the thermodynamic limit. Nevertheless, 
DQMC has the advantage of incorporating non-local cor- 
relations which act in favor of gap formation and DMFT 
overestimates the value of U c because of lacking these 
correlations. 

Due to the single site nature of the DMFT, we were 
not able to address the physics of magnetically ordered 
phases for triangular lattices as one approaches the T = 
limit. These studies have already been performed us- 
ing the DQMC techniquesii2*A^ However, the combina- 
tion of the dynamical cluster approximation (DCA) and 
Qjy[Q30 J 3i 4 32 can a i g0 ^g considered a promising future 
candidate for these studies as a complementary approach 
to DQMC. Work on different magnetically ordered phases 
at different coupling regimes, includin g fer romagnetism 
away from half-filling as discussed in Refl2fjis in progress. 
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